// STL
#include <cmath>                // std::sqrt()
#include <utility>              // std::pair<>, ::make_pair()

// stats++
#include "statsxx/sampling.hpp" // statsxx::sampling::sample_size_m()


//
// DESC: Determine the uncertainty for simultaneous estimation of the parameters of a multinomial population within a distance of the true values at a sample size.
//
// -----
//
// INPUT:
//
//     alpha : significance level
//     n     : sample size
//
// -----
//
// OUTPUT:
//
//     u     : uncertainty
//     m     : the number of nonzero parameters, at which the "worst case" occurs for the given alpha
//
// -----
//
// NOTE: (?!I think!?) The approach *should* not depend on alpha.
//
// NOTE: The method of S. K. Thompson is formulated (only) in terms of alpha, however.
//
// NOTE: ... This is necessary (only) for calculation of m.
//
// -----
//
// NOTE: Using the normal approximation, the success probability is estimated as
//
//     p +/- z*u
//
//     where p is the mean, z is a quantile of a standard normal distribution, and u is the uncertainty.
//
// NOTE: For the binomial distribution, see:
//
//     https://en.wikipedia.org/wiki/Binomial_proportion_confidence_interval
//
// -----
//
// NOTE: See:
//
//     S. K. Thompson, "Sample Size for Estimating Multinomial Proportions," The American Statistician 41, 42--46 (1987)
//
// -----
//
inline std::pair<
                 double, // u
                 int     // m
                 > statsxx::sampling::uncertainty_multinomial(
                                                              const double alpha,
                                                              // -----
                                                              const int    n
                                                              )
{
    double u;
    int    m = statsxx::sampling::sample_size_m(
                                                alpha
                                                );

    // NOTE: There are two ways to calculate this (both described below).
    //
    // -----
    //
    // NOTE: From the S. K. Thompson approach, calculated is:
    //
    //     p +/- d
    //
    // ... (at a given level of alpha.)
    //
    // NOTE: ... We therefore need to calculate:
    //
    //     u = d/z
    //
    // ... (see above).
    //
    // -----
    //
    // METHOD 1:
    //
    // (i) Calculate d2n at the level alpha:
    //
    //     statsxx::sampling::sample_size_multinomial()
    //
    // (ii) Calculate d (using n).
    //
    // (iii) Calculate z:
    //
    //     boost::math::normal dist
    //     z = boost::math::quantile( dist, (1. - alpha/(2*m)) )
    //
    // (iii) Calculate u.
    //
    // -----
    //
    // METHOD 2:
    //
    // (i) Calculate d/z directly.
    //
    // ... Eq. (1) in S. K. Thompson.
    //
    // -----
    //
    // NOTE: METHOD 2 is computationally more straightforward.
    //

    u = std::sqrt( ((1./m)*(1. - 1./m)/n) );

    // RETURN
    return std::make_pair(
                          u,
                          m
                          );
}

//
// ADDITIONAL INFORMATION
//
// =====
//
// NOTE: For binomial data (proportions):
//
// NOTE: ... alpha describes the probability of observing a value outside of (1 - alpha)% of the data.
//
// NOTE: ... With the normal approximation, half of this area is to the left, and the other half to the right.
//
// NOTE: ... Therefore, when calculating the quantile, we would take (1 - alpha/2).
//
// NOTE: For multinomial data (according to Thompson):
//
// NOTE: ... This will be reduced by a factor of m, (1 - alpha/2m) (in the "worst case" scenario).
//
// -----
//
// NOTE: To convert to the standard error:
//
//     SE = (upper limit - lower limit)/(2*z_alpha)
//
//        = (mu + d - (mu - d))/(2*z_alpha)
//
//        = 2*d/(2*z_alpha)
//
//        = d/z_alpha
//
